{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "X_train_path = 'X_train.csv'\n",
    "\n",
    "X = pd.read_csv(X_train_path, index_col=\"id\")\n",
    "\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "scaler = MinMaxScaler(feature_range=(-100, 100))\n",
    "\n",
    "# Scale each row individually for X_train\n",
    "X = X.apply(lambda row: pd.Series(scaler.fit_transform(row.values.reshape(-1, 1)).ravel()), axis=1)\n",
    "\n",
    "# Scale each row individually for X_test\n",
    "X.shape\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import biosppy.signals.ecg as ecg\n",
    "import numpy as np\n",
    "\n",
    "def extract_ecg_features(data, i):\n",
    "    row = data.loc[i].dropna().to_numpy(dtype='float32')\n",
    "    ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate = ecg.ecg(signal=row, sampling_rate=300, show=False)\n",
    "\n",
    "    return ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate\n",
    "\n",
    "\n",
    "ecg_features = ['ts', 'filtered', 'rpeaks', 'templates_ts', 'mean_templates', 'heart_rate_ts', 'heart_rate']\n",
    "ecg_df = pd.DataFrame(columns=ecg_features)\n",
    "ecg_df.index.name = 'id'\n",
    "for i in range(X.shape[0]):\n",
    "    ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate = extract_ecg_features(X, i)\n",
    "    mean_templates = np.mean(templates, axis=0)\n",
    "    ecg_df.loc[len(ecg_df)] = [ts, filtered, rpeaks, templates_ts, mean_templates, heart_rate_ts, heart_rate]\n",
    "    print(i, flush=True)\n",
    "\n",
    "ecg_df.to_csv(\"features/ecg_data.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.fft import fft\n",
    "\n",
    "def extract_fft(data, fs=300):\n",
    "     # Compute the FFT of the signal\n",
    "    if(data.shape[0] == 0):\n",
    "        return np.array([])\n",
    "    fft_values = fft(data)\n",
    "    fft_values = 2.0*np.abs(fft_values[:fs//2])/len(data)\n",
    "\n",
    "    return fft_values\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "fft_features = ['fft_' + e for e in ecg_features]\n",
    "fft_df = pd.DataFrame(columns=fft_features, index=ecg_df.index)\n",
    "for ec in ecg_features:\n",
    "    fft_df['fft_' + ec] = ecg_df[ec].map(extract_fft)\n",
    "    print(ec, end='\\r')\n",
    "\n",
    "fft_df.to_csv(\"features/fft_data.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import welch\n",
    "\n",
    "def extract_psd(data, fs=300):\n",
    "    if(data.shape[0] == 0):\n",
    "        return np.array([])\n",
    "    freqs, psd_values = welch(data, fs)\n",
    "    return psd_values\n",
    "\n",
    "psd_features = ['psd_' + e for e in ecg_features]\n",
    "psd_df = pd.DataFrame(columns=psd_features, index=ecg_df.index)\n",
    "for ec in ecg_features:\n",
    "    psd_df['psd_' + ec] = ecg_df[ec].map(extract_psd)\n",
    "    print(ec, end='\\r')\n",
    "\n",
    "psd_df.to_csv(\"features/psd_data.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_data = pd.concat([ecg_df, fft_df, psd_df], axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import skew, kurtosis\n",
    "\n",
    "def extract_time_features(data, signal_name):\n",
    "    features = {}\n",
    "    if(type(data) != np.ndarray):\n",
    "        return features\n",
    "    if(data.shape[0] == 0):\n",
    "        return features\n",
    "\n",
    "    # Initialize dictionary to hold features\n",
    "\n",
    "    # Basic statistics\n",
    "    features[signal_name + '_mean'] = np.mean(data)\n",
    "    features[signal_name + '_median'] = np.median(data)\n",
    "    features[signal_name + '_std'] = np.std(data)\n",
    "    features[signal_name + '_var'] = np.var(data)\n",
    "    features[signal_name + '_max'] = np.max(data)\n",
    "    features[signal_name + '_min'] = np.min(data)\n",
    "    features[signal_name + '_rms'] = np.sqrt(np.mean(np.square(data)))\n",
    "    features[signal_name + '_peak_to_peak'] = features[signal_name + '_max'] - features[signal_name + '_min']\n",
    "    features[signal_name + '_skewness'] = skew(data)\n",
    "    features[signal_name + '_kurtosis'] = kurtosis(data)\n",
    "\n",
    "    # Zero Crossing Rate\n",
    "    features[signal_name + '_zero_crossing_rate'] = ((data[:-1] * data[1:]) < 0).sum()\n",
    "\n",
    "    # Signal Magnitude Area\n",
    "    features[signal_name + '_sma'] = np.sum(np.abs(data))\n",
    "\n",
    "    # Energy\n",
    "    features[signal_name + '_energy'] = np.sum(np.square(data))\n",
    "\n",
    "    # Entropy\n",
    "    p_signal = data / np.sum(data)  # normalize signal\n",
    "    features[signal_name + '_entropy'] = -np.sum(p_signal*np.log2(p_signal))\n",
    "\n",
    "    # Crest Factor\n",
    "    features[signal_name + '_crest_factor'] = features[signal_name + '_max'] / features[signal_name + '_rms']\n",
    "\n",
    "    # Impulse Factor\n",
    "    features[signal_name + '_impulse_factor'] = features[signal_name + '_max'] / features[signal_name + '_mean']\n",
    "\n",
    "    # Shape Factor\n",
    "    features[signal_name + '_shape_factor'] = features[signal_name + '_rms'] / (np.sum(np.abs(data)) / len(data))\n",
    "\n",
    "    # Clearance Factor\n",
    "    features[signal_name + '_clearance_factor'] = features[signal_name + '_max'] / np.sqrt(np.mean(np.square(np.abs(data))))\n",
    "\n",
    "    # Return the extracted features\n",
    "    return features\n",
    "\n",
    "\n",
    "time_df = pd.DataFrame(index=ecg_df.index)\n",
    "for ec in all_data.columns.to_list():\n",
    "    tf_features = all_data[ec].apply(lambda x : extract_time_features(x, ec)).apply(pd.Series)\n",
    "    time_df = pd.concat([time_df, tf_features], axis=1)\n",
    "    print(ec, end='\\r')\n",
    "        \n",
    "time_df.to_csv(\"features/time_features.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gmean\n",
    "from scipy.stats import entropy\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "def extract_frequency_features(data, signal_name, fs=300):\n",
    "    features = {}\n",
    "    if(type(data) != np.ndarray):\n",
    "        return features\n",
    "    if(data.shape[0] == 0):\n",
    "        return features\n",
    "\n",
    "    # Compute the power spectral density of the signal\n",
    "    freqs, psd_values = welch(data, fs=fs)\n",
    "\n",
    "    # Compute the peak frequency\n",
    "    features[signal_name + '_peak_freq'] = freqs[np.argmax(psd_values)]\n",
    "\n",
    "    # Compute the median frequency\n",
    "    features[signal_name + '_median_freq'] = freqs[len(freqs)//2]\n",
    "\n",
    "    # Compute the spectral entropy\n",
    "    features[signal_name + '_spectral_entropy'] = entropy(psd_values)\n",
    "\n",
    "    # Compute the spectral energy\n",
    "    spectral_energy = np.sum(psd_values)\n",
    "    features[signal_name + '_spectral_energy'] = spectral_energy\n",
    "\n",
    "    # Compute the spectral centroid\n",
    "    spectral_centroid = np.sum(freqs * psd_values) / np.sum(psd_values)\n",
    "    features[signal_name + '_spectral_centroid'] = spectral_centroid\n",
    "\n",
    "    # Compute the spectral spread\n",
    "    spectral_spread = np.sqrt(np.sum((freqs - spectral_centroid)**2 * psd_values) / np.sum(psd_values))\n",
    "    features[signal_name + '_spectral_spread'] = spectral_spread\n",
    "\n",
    "    # Compute the spectral flatness\n",
    "    features[signal_name + '_spectral_flatness'] = gmean(psd_values) / np.mean(psd_values)\n",
    "\n",
    "\n",
    "    # Return the extracted features\n",
    "    return features\n",
    "\n",
    "frequency_df = pd.DataFrame(index=ecg_df.index)\n",
    "\n",
    "for ec in all_data.columns.to_list():\n",
    "    ff = all_data[ec].apply(lambda x : extract_frequency_features(x, ec)).apply(pd.Series)\n",
    "    frequency_df = pd.concat([frequency_df, ff], axis=1)\n",
    "    print(ec, end='\\r')\n",
    "\n",
    "frequency_df.to_csv(\"features/frequency_features.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_poincare_descriptors(data, signal_name):\n",
    "    features = {}\n",
    "\n",
    "    if(type(data) != np.ndarray):\n",
    "        return features\n",
    "    if(data.shape[0] == 0):\n",
    "        return features\n",
    "\n",
    "    # Calculate the differences between successive RR intervals\n",
    "    rr_diff = np.diff(data)\n",
    "\n",
    "    # Calculate SD1 and SD2\n",
    "    features[signal_name + '_SD1'] = np.sqrt(np.std(rr_diff) ** 2 * 0.5)\n",
    "    features[signal_name + '_SD2'] = np.sqrt(2 * np.std(data) ** 2 - 0.5 * np.std(rr_diff) ** 2)\n",
    "\n",
    "    # Calculate SD1/SD2 ratio\n",
    "    features[signal_name + '_SD1/SD2'] = features[signal_name + '_SD1'] / features[signal_name + '_SD2']\n",
    "\n",
    "    # Calculate the area of the Poincaré ellipse\n",
    "    features[signal_name + '_ellipse_area'] = np.pi * features[signal_name + '_SD1'] * features[signal_name + '_SD2']\n",
    "\n",
    "    return features\n",
    "\n",
    "poincarre_df = pd.DataFrame(index=ecg_df.index)\n",
    "\n",
    "for ec in all_data.columns.to_list():\n",
    "    ff = all_data[ec].apply(lambda x : calculate_poincare_descriptors(x, ec)).apply(pd.Series)\n",
    "    poincarre_df = pd.concat([poincarre_df, ff], axis=1)\n",
    "    print(ec, end='\\r')\n",
    "poincarre_df.to_csv(\"features/poincarre_features.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pywt\n",
    "import pandas as pd\n",
    "\n",
    "def extract_wavelet_features(signal, signal_name, wavelet):\n",
    "    features = {}\n",
    "    if type(signal) != np.ndarray:\n",
    "        return features\n",
    "    if signal.shape[0] == 0:\n",
    "        return features\n",
    "\n",
    "    # Compute wavelet coefficients\n",
    "    coeffs = pywt.wavedec(signal, wavelet, level=3, mode='periodic')\n",
    "\n",
    "    # Extract features from each level of coefficients\n",
    "    for i, coeff in enumerate(coeffs):\n",
    "        features.update(extract_time_features(coeff, f\"{signal_name}_{wavelet}_coeff_{i}\"))\n",
    "\n",
    "    return features\n",
    "\n",
    "# Assuming 'signals' is a list of your signals and 'signal_names' is a list of corresponding names\n",
    "wavelets = ['db4']\n",
    "wavelets_df = pd.DataFrame(index=ecg_df.index)\n",
    "\n",
    "\n",
    "for signal_name in all_data.columns.to_list():\n",
    "    for wavelet in wavelets:\n",
    "        features = all_data[signal_name].map(lambda s : extract_wavelet_features(s, signal_name, wavelet)).apply(pd.Series)\n",
    "        wavelets_df = pd.concat([wavelets_df, features], axis=1)\n",
    "        print(signal_name, wavelet, end='\\r')\n",
    "        \n",
    "# Write the DataFrame to a CSV file\n",
    "wavelets_df.to_csv('features/wavelet_features.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_morphological_features(data, signal_name):\n",
    "    features = {}\n",
    "    if(type(data) != np.ndarray):\n",
    "        return features\n",
    "    if (data.shape[0] == 0):\n",
    "        return features\n",
    "\n",
    "    # Compute the R-peak amplitude\n",
    "    r_amplitude = np.max(data)\n",
    "    features[signal_name + '_r_amplitude'] = r_amplitude\n",
    "    r_index = np.argmax(data)\n",
    "    features[signal_name + '_r_index'] = r_index\n",
    "\n",
    "    # Compute the Q and S wave amplitudes\n",
    "    if np.argmax(data) != 0:\n",
    "        q_amplitude = np.min(data[:np.argmax(data)])\n",
    "        features[signal_name + '_q_amplitude'] = q_amplitude\n",
    "        q_index = np.argmin(data[:np.argmax(data)])\n",
    "        features[signal_name + '_q_index'] = q_index\n",
    "\n",
    "    if np.argmax(data) < data.shape[0]:\n",
    "        s_amplitude = np.min(data[np.argmax(data):])\n",
    "        features[signal_name + '_s_amplitude'] = s_amplitude\n",
    "        s_index = np.argmin(data[np.argmax(data):])\n",
    "        features[signal_name + '_s_index'] = s_index\n",
    "\n",
    "    # Compute the T wave amplitude\n",
    "    \n",
    "    if np.argmax(data) < data.shape[0]:\n",
    "        t_amplitude = np.max(data[np.argmax(data):])\n",
    "        features[signal_name + '_t_amplitude'] = t_amplitude\n",
    "        t_index = np.argmax(data[np.argmax(data):])\n",
    "        features[signal_name + '_t_index'] = t_index\n",
    "\n",
    "    # Compute the P wave amplitude\n",
    "    if np.argmin(data) > 0:\n",
    "        p_amplitude = np.max(data[:np.argmin(data)])\n",
    "        features[signal_name + '_p_amplitude'] = p_amplitude\n",
    "        p_index = np.argmax(data[:np.argmin(data)])\n",
    "        features[signal_name + '_p_index'] = p_index\n",
    "\n",
    "    # Compute the QRS complex duration\n",
    "    if np.argmax(data) > 0:\n",
    "        qrs_duration = s_index - q_index\n",
    "        features[signal_name + '_qrs_duration'] = qrs_duration\n",
    "\n",
    "    # Compute the QT interval\n",
    "    if np.argmax(data) < data.shape[0] and np.argmax(data) > 0:\n",
    "        qt_interval = np.argmax(data[np.argmax(data):]) - np.argmin(data[:np.argmax(data)])\n",
    "        features[signal_name + '_qt_interval'] = qt_interval\n",
    "\n",
    "    # Compute the PR interval\n",
    "    if(np.argmin(data) > 0):\n",
    "        pr_interval = r_index - p_index\n",
    "        features[signal_name + '_pr_interval'] = pr_interval\n",
    "\n",
    "    # Compute the ST segment\n",
    "    if np.argmax(data) < data.shape[0] and np.argmax(data) > 0:\n",
    "        st_segment = t_index - s_index\n",
    "        features[signal_name + '_st_segment'] = st_segment\n",
    "\n",
    "    # Compute the P wave duration\n",
    "    if np.argmin(data) > 0:\n",
    "        p_duration = np.argmax(data) - np.argmax(data[:np.argmin(data)])\n",
    "        features[signal_name + '_p_duration'] = p_duration\n",
    "\n",
    "    # Compute the T wave duration\n",
    "    if np.argmax(data) < data.shape[0] and np.argmax(data) > 0:\n",
    "        t_duration = np.argmax(data[np.argmax(data):]) - np.argmax(data)\n",
    "        features[signal_name + '_t_duration'] = t_duration\n",
    "\n",
    "    # Return the extracted features\n",
    "    return features\n",
    "\n",
    "\n",
    "morphological_df = pd.DataFrame(index=ecg_df.index)\n",
    "\n",
    "for ec in all_data.columns.to_list():\n",
    "    mf = all_data[ec].apply(lambda x : extract_morphological_features(x, ec)).apply(pd.Series)\n",
    "    morphological_df = pd.concat([morphological_df, mf], axis=1)\n",
    "    print(ec, end='\\r')\n",
    "morphological_df.to_csv(\"features/morphological_features.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "time_df = pd.read_csv(\"features/time_features.csv\", index_col=\"id\")\n",
    "frequency_df = pd.read_csv(\"features/frequency_features.csv\", index_col=\"id\")\n",
    "poincarre_df = pd.read_csv(\"features/poincarre_features.csv\", index_col=\"id\")\n",
    "wavelets_df = pd.read_csv(\"features/wavelet_features.csv\", index_col=\"id\")\n",
    "morphological_df = pd.read_csv(\"features/morphological_features.csv\", index_col=\"id\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "X_features = pd.concat([time_df, frequency_df, poincarre_df, morphological_df, wavelets_df], axis=1)\n",
    "print(X_features.shape)\n",
    "X_features = X_features.replace([np.inf, -np.inf], np.nan)\n",
    "\n",
    "X_features = X_features.loc[:, ((X_features != 0) & (X_features.notna())).any(axis=0)]\n",
    "print(X_features.shape)\n",
    "\n",
    "# Save the normalized dataframe to a CSV file\n",
    "X_features.to_csv(\"features/X_train_features.csv\")\n",
    "\n",
    "print(X_features.shape)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
